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Abstract. - The dynamics of the 2D Coulomb glass model is investigated by kinetic Monte 
Carlo simulation. An exponential divergence of the relaxation time signals a zero-temperature 
freezing transition. At low temperatures the dynamics of the system is glassy. The local charge 
correlations and the response to perturbations of the local potential show aging. The dynamics 
of formation of the Coulomb gap is slow and the density of states at the Fermi level decays 
in time as a power law. The relevance of these findings for recent transport experiments in 
Anderson-insulating films is pointed out. 



A low-temperature glassy phase that had been anticipated by several authors [1,2] was 
recently observed in studies of hopping conduction in disordered InO films [3,4]. The experi- 
ments show sluggish, non-exponential, relaxation of the conductance [3] as well as aging and 
memory effects [4] similar to those observed in glasses [5]. The details of the dependence of 
the observed effects on carrier density and sample conductance suggests that the glassiness of 
these systems is associated to the interplay of disorder and electron interactions [4]. 

The experimental systems, with very low carrier densities and very large resistances, are 
deep in the Anderson-Mott localized phase. In this regime quantum effects are qualitatively 
unimportant and a description based on the classical 2D Coulomb glass model is appropriate. 
The Hamiltonian of the model is [10] 

Here, r*j denotes the position of the i-th localized state in the plane of the sample, tpi its energy 
e the electron charge and k the background dielectric constant. The electron occupancy 
rii = 0, 1, and there is a uniform compensating background charge density K = l/NJ2i n i 
that ensures global charge neutrality. The presence of the long range Coulomb interaction in 

is a salient feature of this model that has strong effects on its low-temperature properties. 

The Coulomb glass model has been intensively investigated by Monte Carlo simulation but 
so far the main emphasis has been on its low-temperature properties at equilibrium [1,6-9]. 
In this Letter we report the results of the first systematic investigation of its off-equilibrium 
dynamics. 
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The time evolution of Anderson insulators takes place through phonon-assisted processes 
in which electrons hop from occupied sites a to empty sites b with a transition rate [10] 



Here, r is a microscopic timescale, I is the spatial extension of the localized wavefunctions, 
6n a = n a — K and AE a {, is the difference in energy between the final and initial states. The 
origin of the exponential factor in Eq. (J2J is the dependence of the electron-phonon tunneling 
matrix element between sites a and b on the distance between them and <£>(x) is the thermal 
factor, $(x) = e~ x for x > 0, $(x) = 1 for x < 0. The positions of the localized states and 
their energies ipi are both random in general. However, we follow previous authors [6, 7] and 
set (fii = thus limiting our discussion to a simplified model with positional disorder only 
Furthermore, we take K = 1/2 endowing the model with particle-hole symmetry, a feature 
that greatly simplifies the analysis of the results [6,7]. The substitution 2m — 1 = a\ = ±1 
makes Hamiltonian QJ formally equivalent to that of a classical disordered antiferromagnet 
with long range 1/R interactions. The dynamics (J2J) would be quite unusual in the magnetic 
context, however ( 2 ). 

We measure distances and energies in units of the average intersite distance ao and the 
Coulomb energy Ec — e 2 / (k<xq) respectively, and we fix I — ao- Using parameters appropriate 
for InO this corresponds to a carrier density «~3x 10 21 cm~ 3 which is within the range of 
values studied experimentally [3,4]. For these values of the parameters Ec ~ 0.1 eV is huge 
compared to the experimental temperatures. 

We performed kinetic Monte Carlo simulations of model Q . The localization centers are 
randomly distributed within a simulation cell of side L containing N e = L 2 /2 electrons. The 
cell is then periodically repeated in 2D space. The contribution to the potential within the 
simulation cell of charges located outside it is computed using the Ewald method. We start the 
simulations from a random electron configuration to mimic a quench from high temperature 
and then let the system evolve at the working temperature T with the dynamics J2J. The 
latter is implemented as follows. An occupied site a is picked at random and a probability 
P (b | a) oc exp (— 2R a i,/£) is assigned to each of the available empty sites b. A destination 
site is then picked from this probability distribution and an attempt is made to move the 
electron on site a to b. The move is accepted with probability $(AE a b/T). A Monte Carlo 
step (MCS) consists of N such attempts. We simulated systems with L —8, 16, 32 and 64 in 
the temperature range 0.01 < T < 0.1. The length of our runs was typically 2 x 10 6 MCS. 
Physical quantities were monitored as a function of time and the results were averaged over 
between 64 and 6400 realizations of the disorder and the initial conditions, depending on size 
and temperature. 

We begin our discussion by considering the local charge correlation function, 



( x )The symmetry breaking field ipi changes the universality class of the problem but the qualitative features 
of the off-equilibrium dynamics are preserved, at least if the random local potential is small enough. See also 
our concluding remarks 

( 2 )Eq. i s a widely used approximate form for the transition rate in hopping-conduction systems. We 
expect other forms for r a _[, to lead to a qualitatively similar dynamics. 




(2) 





(4) 
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where the brackets denote the average over configurational disorder and thermal noise and 
the waiting time t w is defined as the time elapsed since the quench. If t w is longer than the 
equilibration time r eq there is time-translation invariance and the correlation function will 
depend on the time difference t only. Otherwise, C is expected to depend on both t and t w . 

The T-dependence of r eq was determined using the method of Bhatt and Young [11]. This 
is based on the comparison of two overlap functions, Q r e P (to) = ( [4/iV Snf(to)Sn l l (to)] 2 ) 
and Qt(to) = ([4/iVYJ j 5nf(2t )Sn°;(to)] 2 }. In these expressions the superscripts a and b 
denote two copies of the system (i.e., two systems with the same realization of the disorder) 
that evolve independently starting from different uncorrelated initial conditions. If to > r eq 
the system is at equilibrium in which case [11] Q r ep(to) = Qt(to). The equilibration time 
may thus be defined as the value of to at which the two curves merge (within the statistical 
uncertainty). Fig. shows the time-dependence of these overlap functions for a system of 
size N — 256 and several temperatures. From a quantitative analysis of these results we find 
[cf. inset to the figure] that r eq (T) ~ exp(To/T) where, restoring units, T ~ Ec/2. The 
exponential divergence of r eq at T = is consistent with a freezing transition at T g = as 
reported by other authors [6,8]. The equilibration time is thus finite for all T ^ but it 
becomes very long for T <C Eq- For T > 0.05 our observation time is longer than r eq but the 
latter far exceeds the simulation time for T < 0.04; in this T-range we can only probe the 
off-equilibrium dynamics. These temperature regimes are illustrated in Fig. ^ that shows the 
auto-correlation function for two temperatures, T = 0.05 and T = 0.02, and several values 
of t w . At the higher temperature the relaxation curves for different values of t w coincide for 
t w > 10 5 MCS ~ T eq indicating that equilibrium was attained in the course of the simulation. 
At the lower temperature, however, the i^-dependence persists throughout the observation 
time and the relaxation becomes increasingly sluggish as the system gets older. Equilibrium 
was obviously not reached within the simulation time in this case. 

The size-dependence of our results in the non-stationary regime is remarkably weak as 
shown in Fig. [21 This is the signature of a very slow growth of the size of ordered domains 
£(i) after a quench, a familiar feature of off-equilibrium glassy systems. An analysis of the 
four-point correlation function C^(Rij) = {SnfSn^njSn^) along the lines of Reference [12] 
gives access to £(t) ( 3 ). As illustrated in the inset to Fig. |3| at low temperature, oc log< 

( 3 )Here we use the simple estimate £ 2 (i) = £\ R^C^ (Rij) / • {Rij) . 
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Fig. 2 Fig. 3 

Fig. 2 - The overlap functions Q rep (open symbols) and Qt (filled symbols). The curves have been 
displaced vertically for the sake of clarity. The temperatures are 0.05, 0.06, 0.07, 0.08 and 0.09, form 
bottom to top. Inset : the T-dependence of the equilibration time in MCS. 

Fig. 3 - Size-dependence of C(t + t w ,t w ) for T = 0.03, t w = 10 5 and L=8,16 and 32. Inset: time- 
dependence of the domain size f determined from the analysis of the four point correlation function 
at T = 0.03. 



reaches only a few lattice spacings in the course of our simulations. 

Let us discuss in more detail the off-equilibrium regime illustrated by the data of Fig. [2(b). 
At equilibrium, the correlation function has the scaling form C(t) ~ C eq (t/T) with In- 
dependent A and r. In the non-stationary regime this expression generalizes to [12, 13] 

C(t + t w ,t w )~t- x C a J^^) . (5) 

The above form for the argument of the scaling function C ag , often used in discussions of 
experimental data in glasses [5], also arises in the exact solution of certain mean- field mod- 
els [14]. Here we adopt the widely used parameterization [5] h(x) — exp [a; 1_ ' 1 /(l — fj,)] . In 
the limit /i — > 1 this expression leads to simple aging, C(t + t w ,t w ) ~ t~ x C ag (t/t w ); for other- 
values of /i it leads to an effective relaxation increasing with age as tfc. We used Eq. JSJ to 
perform a scaling analysis of the data of Fig. 2] (b) as well as those obtained for two other 
temperatures, T = 0.01 and T = 0.03. The results are displayed in Fig. 01 We see that Eq. (JSJ 
leads to an excellent collapse of the scaled data for the three temperatures. The exponent A, 
small in magnitude, varies substantially with temperature. The values obtained are consistent 
with limT^o A(T) = as expected for system that freezes at T = 0. The variation of the 
parameter fi with T is much less pronounced. We find simple aging (/i — > 1) as T — + and 
a slight tendency to sub-aging with increasing T. It is worth noticing that simple aging was 
observed in the conductivity of Anderson- insulating films [4] . 

It has been recently suggested [15] that some of the experimentally observed features may 
be due to a slow response of the system to the changes in the local random potential tpi 
[cf. Eq. Q]. We thus studied the response of the system to random perturbations of the 
form Sifi = eiifo where (po 1 is the overall scale of the perturbation and e{ are normalized 
random variables uncorrelated from site to site, {titj) = Sij. This perturbation is switched on 
at time t w and its effect on the system is observed a time t later. The quantity conjugated to 
the random potential is Sn(t) = 1/N £\ (Srii(t)ei). The results of our measurements of the 
response of the system are shown in Fig. JSJ for the case of a random potential of amplitude 



D. R. Grempel: Off-equilibrium dynamics of the two-dimensional Coulomb glass 5 



* 0.8- 

■U 
+ 

■u 

O 0.6 - 

■u 

0.4 L 
10 

Fig. 4 



T 


X 


H 


0.01 


0.003 


~1 


0.02 


0.007 


0.95 


0.03 


0.15 


0.91 



\ 



\ 



10" 10 10 

[<t+t ) ^ - t ^]/(l-n) 



•u 
+ 

•u 




0.1 ■ 



0,05 - 



10 



10 



r=o . 02 

(p =0.01 



10 



10 



Fig. 5 



10 



o t w -0 

□ 10 

o io 2 

A io 2 

v 10 S 

< io s 

> io 6 

10 s 



Fig. 4 - Scaling plot of C(t+t w , t w ). The symbols represent the data obtained for three temperatures, 
T — 0.01,0.02 and 0.03. Inset: the T-dependence of the exponents fj, and A. 



Fig. 5 - Response to a random local field of amplitude (po = 0.01 for T = 0.02. 



(fa = 0.01 and T — 0.02. It may be seen that the response Sn exhibits aging just as the 
correlation function C does. As t w increases, the response gets more and more sluggish and, 
within the time window explored, Sn(t) ~ fni for t^$>t w . 
It can be easily shown that, in the linear response regime, 

Sn{t) = ip Q X (t + t w ,t w ) , X (t, t') = J dt"R(t, t") , (6) 

where x is the local susceptibility and R is the local response function. At equilibrium, x an d 
C are related by the fluctuation dissipation theorem (FDT), T\{t) = 1 — C(t). We show in 
Fig. H3 (a) a parametric plot of the product T% as a function of C for T = 0.02 and several 
waiting times. It can be seen that our data violate this relationship: whereas for each t w the 
points do align on the FDT straight line for short times (C ~ 1), they deviate from it at long 
times (small C). The value of C at which deviations first appear slowly increases with the 
age of the system. 

A generalization of FDT was found to hold for many model systems [14] in which, for long 
t and t', the off-equilibrium x{t,t') still depends on the times only through its dependence 
on C, x{t>t') = 6[C(M')1 where 9 is in general non linear. In glassy systems below the 
glass transition the FDT curves for different values of t w converge to the limiting O [C] as 
t w — > oo. Our data are compatible with [C] consisting of two straight lines with different 
slopes with a break point at some value C* of the correlation. This type of behavior behavior 
is found in some mean-field spin-glass models [14] and in models of structural glasses [17] but 
not in short-range spin glasses. Since our two-dimensional model has no finite-temperature 
glass transition, we should see the break point continuously shift to lower values of C with 
increasing t w such that for t w > r eq the classical FDT line is recovered . This process is slow, 
however, and our data only show a hint of it: on the time scale of the simulation the system 
behaves as a glass that violates FDT. We display in Fig. 0(b) the temperature-dependence of 
the FDT diagram. It is seen that the break point gradually shifts to longer time scales with 
increasing T until, by the time we reach T — 0.08, FDT is fully restored. 

In the study of models of glasses it has proven conceptually useful to introduce an effec- 
tive temperature [16] 1/T cff = —d<d/dC that depends on timescale and represents (loosely 
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Fig. 6 - (a) Parametric plot of susceptibility vs. correlation for T = 0.02 and the waiting times 
indicated in the legend, (b) % vs - C f° r tw = 10 5 and the values of T given in the legend. The 
straight lines have slopes l/T. FDT holds for the highest temperature. T g ~ 0.08 in this T-range. 



speaking) the temperature of the un-equilibrated degrees of freedom. The data displayed in 
Fig- EH (b) suggest that at low T, in the time-domain explored, there exist just two timescales, 
with effective temperatures T e ff = T and T e ff ~ 0.08, respectively. 

Slow relaxation can also be seen in other physical quantities such as the density of states 
p(e) = ^/N(J2i S(e — £;)}. It was suggested [18] that the slow relaxation of the conductance 
observed experimentally may reflect a slow formation of the Efros-Shklovskii Coulomb gap [10] . 
We display in Fig. J7J our results for p(e) for a system of size N — 2500 at T = 0.02. It is 
seen that the density of states, featureless right after the quench, slowly develops a pseudo 
gap for |e| — > 0. The inset to the figure shows the time dependence of the density of states at 
the Fermi level, p(Q, t) ~ t~ v with v w 0.3. A power-law decay of p(0, t) was predicted by the 
phenomenological approach of Reference [18]. 

The glassy behavior described above only exists on timescales shorter than r eq . At fixed 
T, this increases exponentially with the Coulomb energy which is determined by the doping, 
Ec oc 71 1 / 3 . For a carrier density n — 10 21 cm~ 3 , Ec ~ 700X which leads to values of r eq that 
are much greater than experimental timescales at the working temperature. The relaxation 
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Fig. 7 - Time evolution of the Coulomb gap for T = 0.02 for the waiting times shown in the legend. 
The size of the system is N = 2500. Inset: power-law decay of the density of states at the Fermi level. 
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time decreases very rapidly with decreasing n implying an absence of glassincss in the more 
lightly doped samples, a trend that was observed in the experiments [4]. It is interesting to 
notice that 2D spin-glass films, which like the 2D Coulomb glasses discussed here do not have 
a finite temperature phase transition, also exhibit aging effects [19]. 

In summary, we studied the dynamics of the 2D Coulomb-glass model at low temperature. 
The equilibration time of the system diverges exponentially for T — ► signaling a transition at 
T g = 0. At finite T, the evolution of the system after a sudden quench from high temperatures 
is similar to that of glasses below T g for times shorter than T oq . In this regime the charge 
autocorrelation function and the response of to local random perturbations show simple aging 
and FDT violations are observed. The density of states close to the Fermi level evolves slowly 
and p(0) has a power-law decay. 

Our results (except those for the response function) were obtained in the absence of the 
symmetry breaking field ipi. We checked explicitly that the aging behavior of C(t,t') persists 
in the presence of a small random local potential (ipo <C 1). We have preliminary evidence 
that the model with ipo ~ 1 also ages, in contradiction with recent claims in the literature [9]. 

I thank Z. Ovadyahu, A. Vaknin, M. Pollak, L. F. Cugliandolo, A. Kolton and D. Dominguez 
for illuminating discussions. This work was supported in part by NSF Grant No. PHY99- 
07949 and by the program ECOS-Sud, project A01E01. 
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